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Abstract 

It is rigorously shown that the fluctuation formula, which is used in simulations 
to calculate the dielectric constant of interaction site models, corresponds to the 
reaction field with an individual site cut-off rather than with the usual molecular 
center of mass truncation. Within the molecular cut-off scheme, a modified reaction 
field is proposed. An infiuence of the truncation effects is discussed and examined 
by actual Monte Carlo simulations for a MCY water model. 
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1 Introduction 



The calculation of dielectric quantities by computer experiment requires an explicit 
consideration of effects associated with the truncation of long-range interactions. The 
concrete success in this direction has been achieved within the reaction field (RF) geometry 
[1-5]. As a result, computer adapted dielectric theories have been proposed [6-10]. In the 
framework of these theories, a bulk dielectric constant can be determined on the basis of 
a fluctuation formula via correlations obtained in simulations for finite samples. However, 
main attention in the previous investigations has been focused on polar systems with the 
point dipole interaction. As is now well established, the model of point dipoles can not 
reproduce adequately features of real polar liquids. 

At the same time, attempts to apply the RF geometry for more realistic interaction 
site (IS) models have also been made [11-13]. However, acting within a semiphenomeno- 
logical approach, it was not understood how to perform the truncation of intermolecular 
potentials. As a consequence, the molecular cut-off and the usual point dipole RF (PDRF) 
have been assumed. Obviously, such an approach includes effects connected with finitcness 
of the molecule inconsistently. Indeed, the interdipolar potential is replaced by site-site 
Coulomb interactions, whereas the RF is remained in its usual form. An additional com- 
plication for IS models consists in a spatial distribution of charges and this fact is not 
taken into account by the standard PDRF geometry. 

In the present paper we propose two alternative approaches to remedy this situation. 
The first one follows from the usual fluctuation formula which is constructed, however, 
on the microscopic operator of polarization density for IS models. This leads to an ISRF 
geometry, where the cut-off radius is applied with respect to individual charges rather 
than to the molecule as a whole. Nevertheless, the molecular cut-off scheme can also 
be acceptable, but the reaction field together with the fiuctuation formula need to be 
corrected. In the second approach a molecular RF (MRF) geometry is proposed and a 
new quadrupole term is identified. On the basis of a MCY water model we show that 
uncertainties of the dielectric quantities can be significant if the standard PDRF geometry 
is used in computer simulations. 
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2 Interaction site reaction field 



We consider an isotropic, classical system of N identical molecules enclosed in volume 
V. The microscopic electrostatic field created by the molecules at point r eV is equal to 

AT a 

Eir) = E E ^aT^^, = [L{r- r')Q{r')dr' , (1) 

i=l a I' ' i I ^ 

where r" denotes the position for charge g„ of ith molecule, Q{r) — J2i,aQa^{'''' ~ ^i) 
the microscopic operator of charge density, L{p) — —V 1/p and the summation extends 
over all molecules and charged sites. For the investigation of dielectric properties, it is 
more convenient to rewrite the electric field (1) in the polarization representation 

E{r) = / T(r - r')P{r')dr' = -^P(r) + lim / T(r - r')P{r')dr' . (2) 
J 3 p->+o J 

V V 

p<\r-r'\ 

Here T(p) = W 1/p is the dipole-dipole tensor, P{r) denotes the microscopic operator 
of polarization density, defined as V'-P(r) — —Q{r), and the singularity limp^oT(p) = 
— 47r/3 S{p)l has been avoided, where I is the unit tensor of the second rank. The both 
charge (1) and polarization (2) representations are equivalent and applicable for infinite 
{N,V ^ oo) systems. 

In simulations, which deal with finite samples, the sum (1) can not be calculated 
exactly taking into account an infinitely large number of terms. Therefore, we must 
restrict ourselves to a finite set of terms in (1) or to a finite range of the integration in (1) 
and (2) for which \r — r'\ < R, where i? is a cut-off radius. Now the following problem 
appears. How to estimate the cut-off field caused by the integration over unaccessible 
region |r — r'| > i?? The solution of this problem has been found for the first time 
for systems with point dipoles in the RF geometry. The result for conducting boundary 
conditions is [7, 8] 

E{r)^El\r)^-fp{r)+lim^ J (T(r - r') + -^)p(rOdr' , (3) 

V, the 
p<\r-r'\<R 

where a cubic finite sample and toroidal boundary conditions (TBC) have been used, so 
that R < /2. The additional term 1/R^ in the right-hand site of (3) describes the 
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RF which is used for an approximation of the real cut-off field. For a pure spherical cut- 

" sc r 

off (SC) without the RF correction, we have Ej^ (r) = / 7(|r — r'\)L{r — r')Q{r')dr' , 

- sc 

where 7(p) = 1 if p < -R and ■j{p) = otherwise. Obviously, that linii^^oo (r) — 

" RF 

]imR^^E^ (r) = E{r). 

Let us perform the spatial Fourier transform T{k) — J dre~^*^'^ •^(^) for arbitrary 
functions Then one obtains 

Ef{k) = L{k)Q{k) , El\k) = -y P(fc) + (T(fc) + 47r^i^l)P(fc) , (4) 
where 

= -47r(l - j,{kR)) ^ , T(fe) = -y (l - s'^) {skk - l) , (5) 

Q{k) — J2i,aQaG~^^'^' — —ik-P{k), k — 27rn/-v^ is one of the allowed wavevectors of 
the reciprocal lattice, n designates a vector with integer components, k — \k\, k — k/k 
and Jo(-^) — sm{z)/z, j^{z) — — cos{z)/z + sin(2;)/2;^ are the spherical Bessel functions of 
zero and first order, respectively. In view of (5), the relations (4) transform into 

Ef{k) = -A^{l-j,{kR))P^{k) , = -4^(1 , (6) 

where Ph{k) = kk-P{k) = \kQ{k)/k'^ is the longitudinal component of the microscopic 
operator of polarization density. 

It is easy to see from (6) that the both functions Ej^ (fe) and Ej^ (fe) tend to the same 
value E{k) — —ATrPi^{k) of the infinite system at i? — > oo {k ^ 0). However, the results 
converge as R^^ for the pure SC scheme, while as R^^ in the RF geometry, i.e., more 
quickly, because a main part of the truncation effects is taken into account by the RF. This 
is very important in our case, where we hope to reproduce features of infinite systems on 
the basis of finite samples. That is why the pure truncation, which is standard for simple 
fiuids with short-range potentials, is generally not recommended for polar systems with 
long-range nature of the dipolar interaction. The influence of the TBC and the difference 
between micro- and canonical ensembles are of order iV~^ ~ R~^ [14] and, therefore, 
they can be excluded from our consideration. It is worth mentioning that electrostatic 
fields are pure longitudinal. They can be defined via the longitudinal component of the 
microscopic operator of polarization density, that is confirmed by Eq. (6). 
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Let us enclose the system in an external electrostatic field -^^(r) . The material relation 
between the macroscopic polarization Phik) = {^-PL(fc)^ in the weak external field and to- 
tal macroscopic field is 47rPL(fc) = [s^{k) — IjEi^^k), where €^{k) denotes the longitudinal 



wavevector-dependent dielectric constant. Applying the first-order perturbation theory 
with respect to yields for rigid molecules VkBTPi^(k) = (^Pj^{k)'Pi^{—k)'j^E^(k), 
where /cb and T are Boltzmann's constant and temperature, respectively, and is the 
equilibrium average in the absence of the external field. Then, taking into account that 
Eh{k) = E^^{k) + (^E^ (fc)^ and eliminating E^i^k), we obtain the fluctuation formula 

s, (k) - 1 9yG, (k) 

e^{k) l + 27yG^{k)j,{kR)/{kR) ^^^1 )■ K ) 



Here G^{k) = (^Pi^{k)-Pi^{—k)j^ j Njji^ is the longitudinal component of the finite-system 
wavevector-dependent Kirkwood factor, y = AnN ij? j^Vk-oT and fj, = = {Y^aQa^il 
denotes the permanent magnitude of molecule's dipole moment. It is necessary to note 
that we consider rigid IS molecules so that effects associated with molecular and elec- 
tronic polarizabilities are not included in our investigation. In the case of it! — > oo, we 
have j^{kR)/{kR) and computer adapted formula (7) reduces to the well-known fiuc- 
tuation formula for macroscopic systems in terms of the infinite-system Kirkwood factor 
g^{k) = liuiR^^G^ik). 

RF 

As was mentioned earlier, the electric field Ej^ in the form (3), (4) as well as the 
fiuctuation formula (7) have been proposed for the first time to investigate polar systems 
of point dipoles [8]. However, acting within a semiphenomenological framework, it was not 
understood how to perform the truncation of the intermolecular potential (/? at attempts 
to extend this formula for IS models. As a result, the molecular cut-off r^j = \ fi — fj\ < R, 
where rj is the center of mass for iih molecule, and the usual PDRF have been suggested 
[11-13]: 

0?.. = > -. ^ — ^ , Tj,- <R . (8) 

a,b I * 3' 

It is essentially to emphasize that the fluctuation formula (7) takes into account flnite- 
ness of the system explicitly by the factor j^{kR)/{kR). As a result, if the system size is 
sufficiently large (terms of order can be neglected), the bulk {N,V — > oo) dielectric 
constant can be reproduced via the finite-system Kirkwood factor G^{k) which depends 



on it! in a characteristic way. However, to achieve this self-consistency in the evaluation 
of the bulk dielectric constant, the equilibrium averaging in G^{k) must be calculated for 
systems with the intermolecular potential which leads exactly to the microscopic electric 

- RF 

field Ej^ (r) (3). As we shall below, the intermolecular potential (8) does not obey this 
condition. 

To derive the exact intermolecular potential in the charge representation, we perform 

R-F f R,F h, 

the inverse Fourier transform (r) = / dkE^ (fc)e and obtain using (6) 

1-^^^^ Jj,{kR)j,{k\r-rmk] ■ (9) 
/ 

Taking into account that ^ j^{kR)j^{kp)(\.k = p/B? ii p < R and is equal to R/ p^ if 
p > R, we have 

< (-) = Y.^ay^, (l - if \r - <| < (10) 

- RF 

and (r) = otherwise, where the first term in the right-hand side is the Coulomb 
field, while the second contribution corresponds to the RF in the IS description. 

In order to understand nature of this field, we consider a spherical cavity of radius R 
with the center at point r, embedded in an infinite conducting medium. Let us place a 
point charge g„ at point in the cavity, so that \r — r\\ < R. The total electric field e"(r) 
at point r consists of the field due to the charge and the field created by induced charges 
located on the surface of the cavity. According to the method of electrostatic images [5] , 
this last field can be presented as the field of an imaginary charge q* — —qaR/\r — r"| 
which is located at point r*" — r — R^{r — r1)/\r — r"|^ outside the sphere. Then 
e^r) = Qair - <)/|r - r^^ + q^r - 0/|r - <|3 = g„(r - rt){l/\r - r'}\^ - 1/R') 
that is completely in line with the term of sum (10). 

- RF 

In the potential representation {Ej^ (r) — — V<?(r)), we obtain ${r) = J2i,a4>iif), 
where (j)f{r) — qa p1 + \p'f' jR? -|- C), p^ — \r — r"| and C is, in general, an arbitrary 
constant which for infinite systems is chosen as ^1\p<x^oo — 0. In our case, according 
to the toroidal boundary conventional, 0"|pa=^ = whence C = —3/2R~^. Then the 
intermolecular potential of interaction is (p.. = J2a,b Q^'Pi i'f'j) = '}la,b(la4>]{K) = Ea,6<^"^ 
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where 

, \r1-r)\>R 

and the site-site cut-off is performed. 

It is easily seen from (11) that the ISRF part \ J2a,b (la%\K — / transforms into 
the usual form —^^-Hj/E? of point dipoles for < R — d only, where d = 2max|5"| 
is the diameter of the molecule and 5" = rf — ri. In the case if the molecular rather 
than the site-site cut-off is applied to the potential (11), this transformation is valid for 
arbitrary Vij < R. Moreover, in the last case the constant C = —3/2 R~^ is canceled owing 
electroneutrality {J2a Qa = 0) of the molecule and wc recover the result (8) of previous work 
[11]. However, the potential of interaction (11) corresponds completely to the conditions 
at which the fluctuation formula (7) is derived. Therefore, this potential, instead of (8), 
must be used in simulations to obtain a correct value for the dielectric constant. 



3 Molecular reaction field 

In the case of point dipoles, where d — > -|-0, —>■ oo provided fi const, both (8) 
and (11) representations are identical and reduced to the well-known result 

^i,- = -M^-T(r,,)./i,.-^ , nj<R (12) 

for the intcrdipolar interaction in the RF geometry. It is easy to see that in the case of IS 
models, the intermolecular potential (8) takes into account effects associated with finite- 
ness of the molecule inconsistently. For example, the intcrdipolar potential is replaced by 
the real site-site Coulomb ones, whereas the reaction field is remained in its usual form 
of point dipoles. From this point of view a natural question of how to improve the RF 
within the molecular cut-off scheme arises. The simplest way to solve this problem lies in 
the following. 

Let us consider the mentioned above spherical cavity, centered now at some fixed point 
r^, in the infinite conducting medium. We place an ith molecule in such a way that all 
sites of the molecule would be located in the cavity. This condition is fulfilled providing 
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— ■''ol < Rd = R — d/2. The potential of a molecular reaction field at point r belonging 
the cavity can be presented, according to the method of electrostatic images, as 



p-p* 



qaR/pt 



E 



Qa 



R^ pf' 



(13) 



where p = r — and = — r^. Differentiating (13) over r at point yields 



dr 



To 



R^ ' 



drdr 



^ro 



avr(r) 



T-0 



R^ ' 



drdrdr 



To 



(14) 



Here /Xj = Y.aQaPl — Y.aQa^l is the dipole moment of ith molecule, which does not 
depend on owing electroneutrality of the molecule, while qf° = J2aQa{^PiPi — 

and gf° are the tensors of quadrupole and octupole moments, correspondingly, of zth 
molecule with respect to r^. The third rank tensor gf" has the following components 

^^J^aQa^^PiaPipPi-r ~ pfiPiJpi + Pi /s^c^j + pI^^^p))- ^ is more convenient to 
present multipoles of higher order with respect to the molecular center of mass. For the 
tensor of quadrupole moment we obtain qf ° = + Wj, where q^ — Y^a Qai^^t^i ~ ^i"^^) 
is the tensor of quadrupole moment of ith molecule with respect to its center of mass, 
Wj = 3(/iijPj+Pj/Xj)— and Pj = r,— r^. It is necessary to underhne that tensor q, is 
split into dynamical = Y,a Qa^i^i conservative J2a Qa^t^^ parts for rigid molecules. 

Putting — Vj and assuming d <^ R, we obtain the energy of jth molecule in the 
MRF of ith molecule 



dr 



1 avr(r 



drdr 



+ 



i?3 



6 R^ 



+ ... , 



(15) 



where multipoles of higher order have been neglected. Finally, using the RF potential 
ffF — i^fF + yields the desired intermolecular potential 



EH 



„ . \rj - 

a,b I * J ' 



i?3 







nj < Rd 



rij > Rd 



(16) 



where equality q:I = has been used. 

The total reaction field, created by all molecules at point r near r^ is 
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where M{Rd) — X^f -^"^ /Uj and Q(-Rd) = X^f -^"^ Qj denote the total dipole and own 
quadrupole moment, respectively, within the sphere of radius Ra and W(it!d) = YlT~^^ Wj. 
In the case of point dipoles, we have Rd R, Qj, gj, . . . — > and the MRF (17) trans- 
forms into M{R)/R^ + W{R)p/R^. This last formula shows that the reaction field of 
finite systems is inhomogencous even for point dipoles. Only for macroscopic {R — > oo) 
systems, we reproduce the well-known homogeneous reaction field M{R)/R? introduced 
by Barker and Watts [3]. For finite IS systems, additional higher multipole terms appear. 
This brings, for example, into existence of the new quadrupole-dipole and quadrupole- 
quadrupole interactions in the intermolecular potential (16). We note that the idea of 
using the higher multipole moments in the RF has been proposed for the first time by 
Friedman [5]. 

However, the modified intermolecular potential (16) still needs to be complemented 
by a self-consistent fluctuation formula as this has already been done in the preceding 
section by the fluctuation formula (7) for the potential of interaction in the site-site cut- 
off scheme (11). Unfortunately, it is not a simple matter to construct fluctuation formulas 
in the molecular cut-off approach. This problem will be considered in further studying. 

The difference in the RF geometry between IS and PD models hes in the distinction 
for their microscopic operators of polarization density. For IS models 

^L(fc) = ^Ee-^*'"''E9Ae-'*^'^^ =^L(fc)-^fcfc:5:u;,e-*=-^^ + ... , (18) 

1^ i=l a ^ i=l 

where Mi^{k) = kY.f=ik •/ZjC is the microscopic operator of polarization density for 
point dipoles and an expansion over small parameter fc-^" has been made [15]. However, 
putting Pi,{k) = Mi^{k) in the microscopic electric field Ej^ (fe) (6) at the very begin- 
ning and taking attempts to perform the inverse Fourier transform, we obtain that the 
corresponding integral is divergent in fc-spacc when k oo. This divergence is involved 
by the specific nature of point dipoles for which the parameter k-51 becomes indetermi- 
nate in the limit k ^ oo because of — > -|-0 and the expansion (18) fails. Therefore, we 
must manipulate with the full operator Phik) to obtain the interdipolar potential (12) 
consequently and let 6°; +0 at the end of the calculation only. 

9 



Since /j, d and q ~ d^, the quadrupole contribution with respect to the dipole term 
is varied in (16) from of order {d/Kf at Vij = to d/i? at Vij = Rd- Therefore, as far 
as the usual intcrmolccular potential (8) is applied in simulations, the dielectric constant 
can not be reproduced with the precision better than ~ d/i?. It is evident that using 
the modified intermolecular potential (16) will lead to the uncertainties of order [d/Ry. 
They decrease at increasing the size of the sample as R""^, i.e., with the same rate as those 
connected with the truncation of the potential. Effects of the octupole and higher order 
multipole contributions into the MRF are of order (d/R)^ and can be ignored. 



4 Applying the ISRF to a MCY water model 

In the previous investigations [11-13], the standard PDRF geometry (8) has been 
applied to actual simulations of the MCY and TIP4P models. As a result, the static, 
frequency-dependent [11, 12] and wavevector-dependent [13] dielectric constant has been 
determined. For these models d = 1.837A and the cut-off radius R = 9.856A has been 
used in the simulations. From the afore said in the preceding section, it is expected that 
the precision of these calculations can not exceed d/i? ~ 20%. We shall show now by 
actual calculations that this prediction indeed takes place. 

As an example we apply the ISRF geometry (11) to the MCY potential [16]. The 
calculations have been performed with the help of Monte Carlo (MC) simulations, details 
of which are similar to those reported earlier [13], at the density of p= 1.0 g/cm^ and at the 
temperature of T = 292 K, i.e., in the same thermodynamic point and yet with the same 
number N — 256 of molecules and cut-off radius R — 9.856A as considered in [11, 13]. 

Our result of the calculation (7) for the longitudinal components of the wavevector- 
dependent infinite-system Kirkwood factor g^{k) and dielectric constant e^{k) obtained 
within the ISRF geometry is presented in Figs. 1 and 2, respectively, as the full circles 
connected by the solid curves. For the purpose of comparison, analogous calculations 
performed previously [13] within the PDRF are also included in these figures (the open 
circles connected by the dashed curves). It is obvious that oscillations observing in the 
shape of g^{k) and e^{k) obtained within the PDRF method are nonphysical and caused 
by the finite molecular size which is assumed to be zero in this approach. At the same 
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time, the ISRF geometry gives the true, more smooth dependencies for the Kirkwood 
factor and dielectric constant because the influence of the finite molecular size is included 
here explicitly. As we can sec from the figures, deviations of values for the wavcvector- 
dependent dielectric quantities obtained using the PDRF from those evaluated within 
the ISRF geometry are significant. These deviations achieve maximal values about 25% 
near k — SA ^ , where the Kirkwood factor has the first maximum. For great wavevector 
values (A; > 6A ^) the both geometries lead to identical results because the infiuence of 
boundary conditions is negligible in this range of k. 

We remark that the wavevector-dependent quantities were calculated directly for the 
discrete set k — n/cmin of grid points accessible in the simulations, where kmin — 0.319A 
and n is an integer number. These quantities are marked in the figures by the symbols. 
To obtain intermediate values between the grid points we have used the cubic spline inter- 
polation for the most smooth dependency, namely, for g-^{k). Then values of B^{k) can be 
evaluated anywhere in the considered domain of /c-space on the basis of the interpolation 
values of gjyk) via Eq. 7. In particular, the first singularity of ej^k) (see Fig. 2a) has 
been investigated in such a way. 

5 Conclusion 

Two alternative methods (ISRF and MRF) to overcome the difficulties associated with 
finiteness of the molecule with respect to the system size have been proposed for IS mod- 
els of polar systems. It has been shown rigorously that the fiuctuation formula, which 
is commonly used for the calculation of the dielectric constant in computer experiment, 
corresponds to the ISRF geometry with the site-site cut-off for Coulomb interaction po- 
tentials. The molecular cut-off scheme leads to the MRF geometry with an additional 
quadrupole term to the well-known PDRF. 

It has been corroborated by actual calculations that the ISRF geometry exhibits to be 
much more efficient with respect to the usual PDRF method for the investigation of the 
dielectric properties of IS models. The modified MRF approach seem to be comparable 
in efficiency with the ISRF geometry. An application of the MRF to practical simulations 
we hope to perform in further studying. 
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Figure captions 

Fig. 1. Longitudinal component of the wavevector-dependent Kirkwood factor for the 
MCY water. The results in the ISRF and PDRF geometries are plotted by the sohd and 
dashed curves, respectively. 

Fig. 2. Longitudinal component of the wavevector-dependent dielectric constant for 
the MCY water. Notations as for fig. 1. The vertical hues indicate positions of a singu- 
larity. 
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